***  PREAMBLE ******************************************************************

	// clear estimates & macros
estimates clear
macro drop _all
	
* set covariates for models
	
	// pooled OLS
global polscov general dementia service shortterm c.spaces##c.spaces##c.spaces procured lov c.logscb_pop##c.scb_meaninc scb_popshare80plus scb_popshare65plus c.shareprivateyr##c.shareprivateyr i.rc_ideology i.year i.mc

	// facility-FE
global fecov general dementia service shortterm c.spaces##c.spaces##c.spaces lov c.logscb_pop##c.scb_meaninc scb_popshare80plus scb_popshare65plus c.shareprivateyr##c.shareprivateyr i.rc_ideology c.y##c.y##c.y

	// pooled OLS, incl. public (identical as polscov, without procured)
global polscovpub general dementia service shortterm c.spaces##c.spaces##c.spaces lov c.logscb_pop##c.scb_meaninc scb_popshare80plus scb_popshare65plus c.shareprivateyr##c.shareprivateyr i.rc_ideology i.year i.mc

	// set font
graph set window fontface "Times New Roman"
 
	// set paths 
global dataset 		"" 	// <- path to + name of dataset
global datapath 	""	// <- path to dataset
global resultspath  	// <- path to where tables & figures are stored (w/o quotation marks)
********************************************************************************



*** FREQUENCIES, COVERAGE, & SHARES ********************************************
use $dataset, clear

* mean coverage:
bys fid: gen focall=_n
qui count
local noobs = r(N)
qui count if focall == 1
local nofac = r(N)

di "observations:" `noobs'
di "facilities:" `nofac'
di "mean coverage:" `noobs' / `nofac' 

* facilities & observations, per operator category
bys fid opcat: gen foc=_n
forvalues n = 0/4 {
qui count if opcat == `n'		// N(opcat = `n')
local obs_`n' = r(N)
qui count if opcat == `n' & foc == 1		// n(opcat = `n')
local fac_`n' = r(N)

di "***********"
di  "opcat `n':"
di "observations:"
di `obs_`n''
di "facilities:"
di `fac_`n''
di "share of all observations:"
di (`obs_`n'' / `noobs' ) * 100
}

* share private by county
foreach c in 1 3 4 5 6 7 8 9 10 12 13 14 17 18 19 20 21 22 23 24 25 {
use $dataset, clear
qui count if county == `c'
local N_all`c' = r(N)
qui count if county == `c' & opcat != 0		//N(private)
local N_priv`c' = r(N)

di "share private in county = `c':"
di `N_priv`c'' / `N_all`c''  	
}


* PROCURED VS. SELF OWNED 
use $dataset, clear
tab procured if opcat != 0
tab procured nonprofit if opcat != 0, column
tab procured opcat if opcat != 0, column


* COUNT No. OF TRANSITIONS 
use $dataset, clear
keep if opcat!=0
bys fid: gen foc=_n
egen sdnonprofit=sd(nonprofit), by(fid)
egen sdforprofitcat=sd(forprofitcat), by(fid)
count if foc == 1 			
local n_all = r(N)	
	// No. facilities transitioning between NP - FP:
count if foc == 1 & sdnonprofit!=0 & sdnonprofit!=.
	local n_sd = r(N)
	// Share facilities transitioning between NP - FP:
di `n_sd' / `n_all'
	// No. facilities transitioning between PL - PE/PT:
count if foc==1 & sdforprofitcat != 0 & sdforprofitcat != .
	local n_sd = r(N)
	// Share facilities transitioning between PL - PE/PT:
di `n_sd' / `n_all'


* number of non-87301
use $dataset, clear
keep if opcat != 0 
bys fid: gen foc=_n
egen maxb87301 = max(b87301), by(fid)
tab b87301
tab maxb87301 if foc == 1
********************************************************************************

	
	
*** Figure 1: Geographic Distribution of Facilities by Provider Ownership 
// note: only nonprofit has legend included
use $dataset, clear
bys mc:  gen fom=_n 

foreach share in sharenonprofit shareregprofit sharepe shareptrade {
  if "`share'" == "sharenonprofit" {
    local legend "symy(*3) symx(*3) size(*2) position (10)) legorder(hilo) legend(	label(2 "None") label(3 "0 %-10 %" ) label(4 "10 %-25 %" ) label(5 "25 %-50 %" ) label(6 "50 %-75 %" ) label(7 "75 %-100 %" ) size(large)"
  }
  else if "`share'" != "sharenonprofit" {
    local legend "off"
  }
		cd  $datapath
	#delimit;
spmap 	`share' using swehexcoord_mun if fom==1, id(mc) 
		clmethod(custom) clbreaks(0 001 10 25 50 75 100) 
		fcolor(gs0*.001 gs0*.2 gs0*.4 gs0*.6 gs0*.8 gs0*1) 
		ocolor(gs8 ...) osize(*.25 ...) ndocolor(gs6) ndsize(thick) ndfcolor(gs6) 
		polygon(data("swehexcoord_reg.dta") fcolor(none) ocolor(black) osize(*2.5)) 
		legend(`legend')
		title("")  name(`share', replace) ;
	# delimit cr
	*graph export "$resultspath/fig1_map`share'.eps", replace	
	}
********************************************************************************



*** Figure 2: Distribution of Service Quality Indicators 

use $dataset, clear
keep if opcat != 0
foreach var in staffratio nurseratio staffedu plan satisfaction {
  if "`var'" != "staffratio" {
    local cwhite "c(white)"
  }
	sum `var'
	scalar m`var'=r(mean)
	local m_`var': di %6.2f scalar(m`var')
histogram `var', xline(`m_`var'', lp(dash) lw(medthick)) bin(50) blc(gs0) bc(gs8) percent ///
	plotregion(style(none) m(l=0 r=2 t=0 b=0)) ylab(0(25)75, labsize(*2.5) grid) ysc(r(0 78)) 	///
	xlab(, labsize(*2.5) grid format(%9.0f)) graphregion(m(l=0 r=5 t=6 b=3)) xsize(2) ysize(2) ///
	xtitle(" ")  ytitle("%", size(*4) `cwhite') name(newhistogram`var', replace)
	*graph export "$resultspath/fig2_histogram`var'.eps", replace	
	}
********************************************************************************


	
*** Table 2: Summary Statistics 
use $dataset, clear
tab rc_ideology, gen(ide)
label variable ide1						"Left rule"				
label variable ide2						"Cross-ideological rule"				
label variable ide3						"Right rule"	

estpost sum np pl pe pt staffratio nurseratio staffedu plan satisfaction general dementia service shortterm spaces procured	scb_pop logscb_pop scb_meaninc scb_popshare80plus scb_popshare65plus ide1 ide2 ide3 lov shareprivateyr if opcat != 0

esttab . using "$resultspath/table2_summary.tex", replace fragment collabels(\multicolumn{1}{l}{{Obs}} \multicolumn{1}{c}{{Mean}} \multicolumn{1}{c}{{Std.Dev.}} \multicolumn{1}{l}{{Min}} \multicolumn{1}{l}{{Max}}) cells("count(fmt(%13.0fc)) mean(fmt(2 2 2 2 1 1 1 1 1 2 2 2 2 1 2 %13.0fc 2 1 2 2 2 2 2 2 1)) sd(fmt(2 2 2 2 1 1 1 1 1 2 2 2 2 1 2 %13.0fc 2 1 2 2 2 2 2 2 1)) min(fmt(0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 %13.0fc 2 1 2 2 0 0 0 0 0)) max(fmt(0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 %13.0fc 2 1 2 2 0 0 0 0 0))") label nonumber noobs  booktabs alignment(S) order(np pl pe pt staffratio nurseratio staffedu plan satisfaction general dementia service shortterm spaces procured scb_pop logscb_pop scb_meaninc scb_popshare80plus scb_popshare65plus ide1 ide2 ide3 lov) refcat(np "\textbf{Provider ownership}" staffratio "\textbf{Outcomes}"  general "\textbf{Facility-Level Controls}" scb_pop "\textbf{Municipal-Level Controls}", nolabel) type
********************************************************************************



*** Figure 3: Provider Ownership and Indicators of Service Quality: Bivariate Relationships ***
use $dataset, clear
foreach var in staffratio nurseratio staffedu plan satisfaction {
  if "`var'" == "staffratio" {
    local yrange "0(10)30"
    local lmargin "l=5"
  }
  else if "`var'" == "nurseratio" {
    local yrange "0(1)6"
    local lmargin "l=10"
  }
  else if "`var'" == "staffedu" | "`var'" == "plan" | "`var'" == "satisfaction" {
    local yrange "50(10)100"
    local lmargin "l=0"
  }
reg `var' i.opcat, cluster(fid)
	margins opcat
	marginsplot, recast(bar) plotop(barw(.5) fc(gs8) lc(gs0)) ciopts(lc(gs0) lw(medthick)) title("") ytitle("") ylab(`yrange',grid labsize(*3) angle(0)) xtitle("") xlab(, labsize(*3) angle(25) labgap(small) notick) xsize(1.5) ysize(2) graphregion(m(`lmargin')) name(`var'bivauto, replace)
*graph export "$resultspath/fig3_`var'biv.eps", replace
}
********************************************************************************

 

 
*** Figure 4: Provider Ownership and Indicators of Service Quality: pooled OLS with municipal and time fixed effects 
use $dataset, clear
drop if opcat == 0
foreach var in staffratio nurseratio staffedu plan satisfaction {
  if "`var'" == "staffratio" {
    local yrange "0(10)30"
    local lmargin "l=5"
  }
  else if "`var'" == "nurseratio" {
    local yrange "0(1)6"
    local lmargin "l=10"
  }
  else if "`var'" == "staffedu" | "`var'" == "plan" | "`var'" == "satisfaction" {
    local yrange "50(10)100"
    local lmargin "l=0"
  }
 reg `var' i.opcat $polscov, cluster(fid)
di "Avg. Adjusted predictions: `var'"
	margins opcat
	marginsplot, recast(bar) plotop(barw(.5) fc(gs8) lc(gs0)) ciopts(lc(gs0) lw(medthick)) title("") ytitle("") ylab(`yrange',grid labsize(*3) angle(0)) xtitle("") xlab(, labsize(*3) angle(25) labgap(small) notick) xsize(1.5) ysize(2) graphregion(m(`lmargin')) name(`var'bivauto, replace)
*graph export "$resultspath/fig4_`var'multi.eps", replace
}
********************************************************************************




*** Table 3: H2: Ownership and Indicators of Service Quality: PL vs. PE and PT, FE estimates
* Producing underlying estimates
use $dataset, clear
drop if opcat == 0

foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
 xtreg `var' i.opcat $fecov, vce(cluster fid) fe
eststo h2fe`var': margins opcat, pwcompare(effects)  post
estadd local fe "$\checkmark$"
estadd matrix beta_vs= r(table_vs)["b", 4..5]
estadd matrix se_vs= r(table_vs)["se", 4..5]
estadd matrix p_vs= r(table_vs)["pvalue", 4..5] 
}
  
* Procuding tables
esttab h2festaffratio h2fenurseratio h2festaffedu h2feplan h2fesatisfaction using "$resultspath/table3_h2fe.tex", replace cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) collabels(none) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") compress nogap  noobs nonumbers stats(fe N, fmt(%13.0fc) labels("Facility FEs"  "Obs.")) booktabs fragment t(2) b(2) 

	// Add no. facilities
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
eststo h2fe`var': qui xtreg `var' i.opcat $polscov, cluster(fid) fe
estadd scalar fac = e(N_clust)  
}
esttab h2festaffratio h2fenurseratio h2festaffedu h2feplan h2fesatisfaction using "$resultspath/table3_h2fe.tex", append  compress nogap drop(*) nomtitles nolines noobs nonumbers stats( fac, fmt(0) labels("Facilities")) booktabs fragment 
********************************************************************************




*** Figure 5: Ownership and Indicators of Service Quality: Alternative Estimation Techniques
use $dataset, clear
tab mc, gen(munn)
tab year, gen(yr)
tab rc_ideology, gen(ide)
gen logscb_popxscb_meaninc = logscb_pop*scb_meaninc
gen shareprivateyr_2 = shareprivateyr*shareprivateyr
gen spaces_2 = spaces*spaces
gen spaces_3 = spaces*spaces*spaces

drop if opcat == 0

* produce estimates
foreach var in staffratio nurseratio staffedu plan satisfaction {
foreach h in 1 2 {
	// Main
reg `var' ib`h'.opcat $polscov, cluster(fid)
	est store h`h'_munfe_`var'
	capture predict `var'sample if e(sample)
	// LDV
reg `var' ib`h'.opcat l.`var' $polscov, cluster(fid)
	est store h`h'_ldv_`var'
	// RE
xtreg `var' ib`h'.opcat $polscov, cluster(fid) 
	est store h`h'_re_`var'
	// BE
xtreg `var' ib`h'.opcat $polscov, be
	est store h`h'_be_`var'
	// HT
xthtaylor `var' pl pe pt general dementia service shortterm spaces spaces_2 spaces_3 procured lov logscb_pop scb_meaninc logscb_popxscb_meaninc ide2 ide3 scb_popshare80plus scb_popshare65plus shareprivateyr shareprivateyr_2  yr* munn*, vce(cluster fid) endog(pl pe pt general dementia service shortterm spaces spaces_2 spaces_3 procured) 
	est store h1_ht_`var'	
xthtaylor `var' np pe pt general dementia service shortterm spaces spaces_2 spaces_3 procured lov logscb_pop scb_meaninc logscb_popxscb_meaninc ide2 ide3 scb_popshare80plus scb_popshare65plus shareprivateyr shareprivateyr_2  yr* munn*, vce(cluster fid) endog(np pe pt general dementia service shortterm spaces spaces_2 spaces_3 procured) 
	est store h2_ht_`var'	
	// FE
xtreg `var' ib`h'.opcat $fecov, vce(cluster fid) fe
	est store h`h'_fe_`var' 
	}
	
* produce coeficient plots	
capture gen str orig="{bf:P}"
capture gen str ldv="{bf:LDV}"
capture gen str re="{bf:RE}"
capture gen str be="{bf:BE}"
capture gen str ht="{bf:HT}"
capture gen str fe="{bf:FE}"

 if "`var'" == "staffratio" {
    local title "Staff density"
 	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "nurseratio" {
    local title "Nurse density"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "staffedu"  {
    local title "Staff education"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "plan"  {
    local title "Updated plan"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "satisfaction" {
    local title "Satisfaction"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }	
  
	// H1
coefplot 	(h1_munfe_`var', mlabel(orig)) 	///
			(h1_ldv_`var'	, mlabel(ldv)) 	///
			(h1_be_`var'	, mlabel(be))	///
			(h1_re_`var'	, mlabel(re)) 	///
	,  mc(gs0) mlabc(white) keep(pl pe pt 2.opcat 3.opcat 4.opcat) rename(pl = 2.opcat pe = 3.opcat pt = 4.opcat) transform(* = @/`sd') xline(0, lw(thick) lc(gs4)) xlab(-1 (.5) .5, labsize(*3) grid) ylab(, nolab notick) xsc(r(-1 .5)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) name(h1_techniques_`var', replace)
	
	// H2
coefplot 	(h2_munfe_`var', mlabel(orig)) 	///
			(h2_ldv_`var'	, mlabel(ldv)) 	///
			(h2_be_`var'	, mlabel(be))	///
			(h2_re_`var'	, mlabel(re)) 	///
			(h2_ht_`var'	, mlabel(ht))	///
			(h2_fe_`var'	, mlabel(fe))	///
	,   mc(gs0) mlabc(white) keep(pe pt 3.opcat 4.opcat) rename(pe = 3.opcat pt = 4.opcat) transform(* = @/`sd') xline(0, lw(thick) lc(gs4)) xlab(-1 (.5) .5, labsize(*3) grid) ylab(, nolab notick) xsc(r(-1 .5)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) name(h2_techniques_`var', replace)
}	

* produce combined figures
	// H1
	capture gen one=1
	capture gen two=2
tw sc one two, plotregion(style(none) m(medium)) m(i) ysc(off) xsc(off) xsize(1) ysize(10) text(1.55 3 "PL", size(*3.5)) text(1 3 "PE", size(*3.5)) text(.45 3 "PT", size(*3.5))  name(h1_technique_title, replace)
gr combine 	h1_technique_title 	h1_techniques_staffratio h1_techniques_nurseratio h1_techniques_staffedu h1_techniques_plan h1_techniques_satisfaction, rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(tech_h1, replace)
*graph export "$resultspath/fig5_tech_h1.eps", replace

	// H2
tw sc one two, plotregion(style(none) m(medium)) m(i) ysc(off) xsc(off) xsize(1) ysize(10) text(1.5 3 "PE", size(*3.5)) text(.5 3 "PT", size(*3.5))  name(h2_techniques_title, replace)
gr combine 	h2_techniques_title h2_techniques_staffratio h2_techniques_nurseratio h2_techniques_staffedu h2_techniques_plan h2_techniques_satisfaction,	rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(tech_h2, replace)
*graph export "$resultspath/fig5_tech_h2.eps", replace
********************************************************************************




*** Figure 6: Ownership and Indicators of Service Quality: Alternate Samples
use $dataset, clear
drop if opcat == 0
* produce estimates
foreach var in staffratio nurseratio staffedu plan satisfaction {
	// original
reg `var' ib1.opcat $polscov, cluster(fid)
	est store h1_orig_`var'
	predict `var'sample if e(sample)
reg `var' ib2.opcat $polscov, cluster(fid)
	est store h2_orig_`var'
	// SNI=87301
reg `var' ib1.opcat $polscov if b87301 == 1, cluster(fid)
	est store h1_87301_`var'
reg `var' ib2.opcat $polscov if b87301 == 1, cluster(fid)
	est store h2_87301_`var'
	// Ekonomisk Förening = PL 
reg `var' ib1.opcatalt $polscov, cluster(fid)
	est store h1_alt_`var'
reg `var' ib2.opcatalt $polscov, cluster(fid)
	est store h2_alt_`var'

* produce coeficient plots	
capture gen str orig="{bf:O}"
capture gen str _87301="{bf:BR}"
capture gen str alt="{bf:EA}"
	
  if "`var'" == "staffratio" {
    local title "Staff density"
 	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "nurseratio" {
    local title "Nurse density"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "staffedu"  {
    local title "Staff education"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "plan"  {
    local title "Updated plan"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "satisfaction" {
    local title "Satisfaction"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }	
  
	// H1
coefplot 	(h1_orig_`var', mlabel(orig)) ///
			(h1_87301_`var', mlabel(_87301)) ///
			(h1_alt_`var', mlabel(alt)), ///
mc(gs0) mlabc(white) keep(2.opcat 3.opcat 4.opcat 2.opcatalt 3.opcatalt 4.opcatalt) rename(2.opcatalt = 2.opcat 3.opcatalt = 3.opcat 4.opcatalt = 4.opcat) transform(* = @/`sd') xline(0, lw(thick) lc(gs4)) xlab(-1 (.5) .5, labsize(*3) grid) ylab(, nolab notick) xsc(r(-1 .5)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) name(h1_samples_`var', replace)

	// H2
coefplot 	(h2_orig_`var', mlabel(orig)) ///
			(h2_87301_`var', mlabel(_87301)) ///
			(h2_alt_`var', mlabel(alt)), ///
mc(gs0) mlabc(white) keep(3.opcat 4.opcat 3.opcatalt 4.opcatalt) rename(3.opcatalt = 3.opcat 4.opcatalt = 4.opcat) transform(* = @/`sd') xline(0, lw(thick) lc(gs4)) xlab(-1 (.5) .5, labsize(*3) grid) ylab(, nolab notick) xsc(r(-1 .5)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) name(h2_samples_`var', replace)
}

* produce combined figures
	// H1
	capture gen one=1
	capture gen two=2
tw sc one two, plotregion(style(none) m(medium)) m(i) ysc(off) xsc(off) xsize(1) ysize(10) text(1.55 3 "PL", size(*3.5)) text(1 3 "PE", size(*3.5)) text(.45 3 "PT", size(*3.5))  name(h1_samples_title, replace)
gr combine 	h1_samples_title 	h1_samples_staffratio h1_samples_nurseratio h1_samples_staffedu h1_samples_plan h1_samples_satisfaction, rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(samples_h1, replace)
*graph export "$resultspath/fig6_sample_h1.eps", replace

	// H2
tw sc one two, plotregion(style(none) m(medium)) m(i) ysc(off) xsc(off) xsize(1) ysize(10) text(1.5 3 "PE", size(*3.5)) text(.5 3 "PT", size(*3.5))  name(h2_samples_title, replace)
gr combine 	h2_samples_title h2_samples_staffratio h2_samples_nurseratio h2_samples_staffedu h2_samples_plan h2_samples_satisfaction,	rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(samples_h2, replace)
*graph export 	"$resultspath/fig6_sample_h2.eps", replace
********************************************************************************






********************************************************************************
************************* ONLINE APPENDIX **************************************
********************************************************************************


*** D. Relationships between quality of care indicators ************************

*** Table D1: Correlations between quality of care indicators
use $dataset, clear
drop if opcat == 0
label variable staffedu "\shortstack{Share of Staff w/ \\\ Appropriate Education}"
local vars staffratio nurseratio staffedu plan satisfaction
pwcorr `vars', star(.05) 
estpost correlate `vars', matrix 
esttab . using "$resultspath/tabled1_dvcorrelations.tex", unstack b(2) p noobs nonumber star(* .05) label booktabs fragment replace type  compress nogap
********************************************************************************




*** E. Comparison of distributions of place-based demographic and socioeconomic confounders in nonprofits and for-profits **************************************

use $dataset, clear
gen logscb_popxscb_meaninc = logscb_pop*scb_meaninc


tw histogram logscb_pop if opcat == 1, color(none) lcolor(red) percent  w(.1)  start(8) || histogram logscb_pop if opcat == 2 | opcat == 3 | opcat == 4, color(none) lcolor(blue) percent w(.1)  start(8) xlab(,grid) ylab(,grid) xtitle("log(Population)") graphregion(m(l=0 r=2.5 t=6 b=6)) xsize(2) ysize(2) legend(order(1 "Nonprofits" 2 "For-profits")) name(logscb_pop, replace)
*graph export "$resultspath/fige1_logscb_pop.eps", replace

tw histogram scb_meaninc if opcat == 1, color(none) lcolor(red) percent  w(8)  start(150) || histogram scb_meaninc if opcat == 2 | opcat == 3 | opcat == 4, color(none) lcolor(blue)  percent w(8) start(150) xlab(,grid) ylab(,grid) xtitle("Mean income") graphregion(m(l=0 r=2.5 t=6 b=6)) xsize(2) ysize(2) legend(order(1 "Nonprofits" 2 "For-profits")) name(scb_meaninc, replace)
*graph export "$resultspath/fige1_scb_meaninc.eps", replace

tw histogram logscb_popxscb_meaninc if opcat == 1, color(none) lcolor(red) percent  w(100)  start(1500) || histogram logscb_popxscb_meaninc if opcat == 2 | opcat == 3 | opcat == 4, color(none) lcolor(blue)  percent w(100) start(1500) xlab(,grid) ylab(,grid) xtitle("Mean income * log(Population)") graphregion(m(l=0 r=2.5 t=6 b=6)) xsize(2) ysize(2) legend(order(1 "Nonprofits" 2 "For-profits")) name(scb_meanincxlogscb_pop, replace)
*graph export "$resultspath/fige1_scb_meanincxlogscb_pop.eps", replace
********************************************************************************




*** F. Results of the analysis omitting facility-level confounders *************

global polscovnofac  lov c.logscb_pop##c.scb_meaninc scb_popshare80plus scb_popshare65plus c.shareprivateyr##c.shareprivateyr i.rc_ideology i.year i.mc
global fecovnofac lov c.logscb_pop##c.scb_meaninc scb_popshare80plus scb_popshare65plus c.shareprivateyr##c.shareprivateyr i.rc_ideology c.y##c.y##c.y

* Producing underlying estimates
	// Pooled estimates
use $dataset, clear
drop if opcat == 0

foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
qui reg `var' i.opcat $polscovnofac, cluster(fid)
eststo h1pols`var': margins opcat, pwcompare(effects) post
estadd matrix beta_vs= r(table_vs)["b", 1..3]
estadd matrix se_vs= r(table_vs)["se", 1..3]
estadd matrix p_vs= r(table_vs)["pvalue", 1..3] 

qui reg `var' i.opcat $polscovnofac, cluster(fid)
eststo h2pols`var': margins opcat, pwcompare(effects)  post
estadd matrix beta_vs= r(table_vs)["b", 4..5]
estadd matrix se_vs= r(table_vs)["se", 4..5]
estadd matrix p_vs= r(table_vs)["pvalue", 4..5] 

	// FE estimates
qui xtreg `var' i.opcat $fecovnofac, vce(cluster fid) fe
eststo h1fe`var': margins opcat, pwcompare(effects) post
estadd matrix beta_vs= r(table_vs)["b", 1..3]
estadd matrix se_vs= r(table_vs)["se", 1..3]
estadd matrix p_vs= r(table_vs)["pvalue", 1..3] 

qui xtreg `var' i.opcat $fecovnofac, vce(cluster fid) fe
eststo h2fe`var': margins opcat, pwcompare(effects)  post
estadd matrix beta_vs= r(table_vs)["b", 4..5]
estadd matrix se_vs= r(table_vs)["se", 4..5]
estadd matrix p_vs= r(table_vs)["pvalue", 4..5] 
}


*** Table F1: H1: Ownership and Indicators of Service Quality: NP vs. for-Profits
esttab h1polsstaffratio h1polsnurseratio h1polsstaffedu h1polsplan h1polssatisfaction using "$resultspath/tablef1_h1nofac.tex", replace cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) collabels(none) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") compress nogap  noobs nonumbers stats(N, fmt(%13.0fc) labels("Obs.")) booktabs fragment t(2) b(2) 

	// Add no. facilities
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
eststo h1pols`var': qui reg `var' i.opcat $polscovnofac, cluster(fid)
estadd scalar fac = e(N_clust)  
}
esttab h1polsstaffratio h1polsnurseratio h1polsstaffedu h1polsplan h1polssatisfaction using "$resultspath/tablef1_h1nofac.tex", append  compress nogap drop(*) nomtitles nolines noobs nonumbers stats(fac, fmt(0) labels("Facilities")) booktabs fragment 



*** Table F2: H2: Ownership and Indicators of Service Quality: PL vs. PE and PT
esttab h2polsstaffratio h2polsnurseratio h2polsstaffedu h2polsplan h2polssatisfaction using "$resultspath/tablef2_h2nofac.tex", replace collabels(none) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") drop(*) compress nogap noobs nonumbers  booktabs fragment

esttab h2polsstaffratio h2polsnurseratio h2polsstaffedu h2polsplan h2polssatisfaction using "$resultspath/tablef2_h2nofac.tex", append cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.10 * 0.05 ** 0.01 *** 0.001) collabels(none) nomtitle eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") mgroups("Panel A: POLS", pattern(1 0 0 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) compress nogap nolines noobs nonumbers  booktabs fragment

esttab h2festaffratio h2fenurseratio h2festaffedu h2feplan h2fesatisfaction using "$resultspath/tablef2_h2nofac.tex", append cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.10 * 0.05 ** 0.01 *** 0.001) collabels(none) nomtitle eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") mgroups("Panel B: FE", pattern(1 0 0 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) compress nogap nolines noobs nonumbers  booktabs fragment stats(N , fmt(%13.0fc) labels("\midrule Obs.")) 

	// Add no. facilities
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
eststo h2pols`var': qui reg `var' i.opcat $polscovnofac, cluster(fid)
estadd scalar fac = e(N_clust)  
}
esttab h2polsstaffratio h2polsnurseratio h2polsstaffedu h2polsplan h2polssatisfaction using "$resultspath/tablef2_h2nofac.tex", append  compress nogap drop(*) nomtitles nolines noobs nonumbers stats(fac, fmt(0) labels( "Facilities")) booktabs fragment 
********************************************************************************




*** G. Results of the analysis using a complete sample, including public facilities
use $dataset, clear

	// count number of observations by opcat that varies from public and private
tab opcat if facsample == 1

*** Figure G1: Results, including publicly operated facilities
* produce regression estimates
	// POLS
foreach var in staffratio nurseratio staffedu plan satisfaction {
reg `var' i.opcat $polscovpub, cluster(fid)	// Full sample
est store polsfull_`var'
capture predict `var'sample if e(sample)
reg `var' i.opcat $polscovpub if munsample == 1, cluster(fid)	// Municipalities with both public & private
est store polsmunsample_`var'
reg `var' i.opcat $polscovpub if facsample == 1, cluster(fid)	// Facilities with both public & private
est store polsfacsample_`var'

	// FE
xtreg `var' i.opcat $fecov, vce(cluster fid) fe	// Full sample
est store fefull_`var'
xtreg `var' i.opcat $fecov if munsample == 1, vce(cluster fid) fe	// Municipalities with both public & private
est store femunsample_`var'
xtreg `var' i.opcat $fecov if facsample == 1, vce(cluster fid) fe	// Facilities with both public & private
est store fefacsample_`var'
}

	// count number of observations per model
	esttab polsfull_* polsmunsample_* polsfacsample_*, drop(*)

* produce figures	
capture gen str full="{bf:F}"
capture gen str pim="{bf:M}"
capture gen str ep="{bf:P}"
foreach var in staffratio nurseratio staffedu plan satisfaction {
  if "`var'" == "staffratio" {
    local title "Staff density"
 	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "nurseratio" {
    local title "Nurse density"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "staffedu"  {
    local title "Staff education"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "plan"  {
    local title "Updated plan"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }
  else if "`var'" == "satisfaction" {
    local title "Satisfaction"
	sum `var' if `var'sample != .
	local sd = r(sd)
	di `sd'
  }	
coefplot 	(polsfull_`var', mlabel(full)) 	///
			(polsmunsample_`var', mlabel(pim)) 	///
			(polsfacsample_`var', mlabel(ep)) 	///
	,  mc(gs0) mlabc(white) keep(1.opcat 2.opcat 3.opcat 4.opcat) xline(0, lw(thick) lc(gs4)) xlab(, labsize(*2) grid) ylab(, nolab notick) xsc(r(0)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) transform(* = @/`sd') name(polspublic_`var', replace)
	
	
coefplot 	(fefull_`var', mlabel(full)) 	///
			(femunsample_`var'	, mlabel(pim)) 	///
			(fefacsample_`var'	, mlabel(ep)) 	///
	,  mc(gs0) mlabc(white) keep(1.opcat 2.opcat 3.opcat 4.opcat) xline(0, lw(thick) lc(gs4)) xlab(, labsize(*2) grid) ylab(, nolab notick) xsc(r(0)) legend(off) ci(95) ciopts(recast(rcap) color(gs0)  lwidth(medium)) m(O) msize(*5) mlabsize(*1.5) mlabpos(0) xsize(6) ysize(10) title(`title', size(*2)) transform(* = @/`sd') name(fepublic_`var', replace)
	}
	
* combined figures

	// POLS
	capture gen one=1
	capture gen two=2
tw sc one two, plotregion(style(none) m(medium)) m(i) ysc(off) xsc(off) xsize(1) ysize(10) text(1.75 3 "NP", size(*3.5)) text(1.25 3 "PL", size(*3.5)) text(.75 3 "PE", size(*3.5)) text(.25 3 "PT", size(*3.5)) name(public_title, replace)

gr combine 	public_title polspublic_staffratio polspublic_nurseratio polspublic_staffedu polspublic_plan polspublic_satisfaction, rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(fullpols, replace)
*graph export 	"$resultspath/figg1_fullpols.eps", replace

	// FE
gr combine 	public_title fepublic_staffratio fepublic_nurseratio fepublic_staffedu fepublic_plan fepublic_satisfaction, rows(1) xsize(37) ysize(10) graphregion(m(t=0 b=0 l=-60 r=0)) name(fullfe, replace)
*graph export "$resultspath/figg1_fullfe.eps", replace
********************************************************************************




*** H. Hypothesis tests, pooled OLS estimates **********************************
* Producing underlying estimates
	// Pooled estimates
use $dataset, clear
drop if opcat == 0

foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
qui reg `var' i.opcat $polscov, cluster(fid)
eststo h1pols`var': margins opcat, pwcompare(effects) post
estadd matrix beta_vs= r(table_vs)["b", 1..3]
estadd matrix se_vs= r(table_vs)["se", 1..3]
estadd matrix p_vs= r(table_vs)["pvalue", 1..3] 

qui reg `var' i.opcat $polscov, cluster(fid)
eststo h2pols`var': margins opcat, pwcompare(effects)  post
estadd matrix beta_vs= r(table_vs)["b", 4..5]
estadd matrix se_vs= r(table_vs)["se", 4..5]
estadd matrix p_vs= r(table_vs)["pvalue", 4..5] 

	// FE estimates
qui xtreg `var' i.opcat $fecov, vce(cluster fid) fe
eststo h1fe`var': margins opcat, pwcompare(effects) post
estadd matrix beta_vs= r(table_vs)["b", 1..3]
estadd matrix se_vs= r(table_vs)["se", 1..3]
estadd matrix p_vs= r(table_vs)["pvalue", 1..3] 

qui xtreg `var' i.opcat $fecov, vce(cluster fid) fe
eststo h2fe`var': margins opcat, pwcompare(effects)  post
estadd matrix beta_vs= r(table_vs)["b", 4..5]
estadd matrix se_vs= r(table_vs)["se", 4..5]
estadd matrix p_vs= r(table_vs)["pvalue", 4..5] 
}
  
*** Table H1: H1: Ownership and Indicators of Service Quality: NP vs. for-profits
esttab h1polsstaffratio h1polsnurseratio h1polsstaffedu h1polsplan h1polssatisfaction using "$resultspath/tableh1_h1.tex", replace cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) collabels(none) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") compress nogap  noobs nonumbers stats(N, fmt(%13.0fc) labels("Obs.")) booktabs fragment t(2) b(2) 

	// Add no. facilities
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
eststo h1pols`var': qui reg `var' i.opcat $polscov, cluster(fid)
estadd scalar fac = e(N_clust)  
}
esttab h1polsstaffratio h1polsnurseratio h1polsstaffedu h1polsplan h1polssatisfaction using "$resultspath/tableh1_h1.tex", append  compress nogap drop(*) nomtitles nolines noobs nonumbers stats(fac, fmt(0) labels("Facilities")) booktabs fragment 


*** Table H2: H2: Ownership and Indicators of Service Quality: PL vs. PE and PT
esttab h2polsstaffratio h2polsnurseratio h2polsstaffedu h2polsplan h2polssatisfaction using "$resultspath/tableh2_h2.tex", replace cells(beta_vs(star pvalue(p_vs) fmt(%9.3f) ) p_vs(par)) starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) collabels(none) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") eqlab(none) rename(3vs2.opcat 3vs1.opcat 4vs2.opcat 4vs1.opcat ) varlabels(2vs1.opcat "PL" 3vs1.opcat "PE" 4vs1.opcat "PT") compress nogap  noobs nonumbers stats(N, fmt(%13.0fc) labels("Obs.")) booktabs fragment t(2) b(2) 

	// Add no. facilities
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
eststo h2pols`var': qui reg `var' i.opcat $polscov, cluster(fid)
estadd scalar fac = e(N_clust)  
}
esttab h2polsstaffratio h2polsnurseratio h2polsstaffedu h2polsplan h2polssatisfaction using "$resultspath/tableh2_h2.tex", append  compress nogap drop(*) nomtitles nolines noobs nonumbers stats(fac, fmt(0) labels("Facilities")) booktabs fragment 
********************************************************************************
	

	
	
*** I. Main results, full ******************************************************

* Producing underlying estimates
use $dataset, clear
drop if opcat == 0
label var y "Year"
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
	// Pooled estimates
eststo baselinepols`var':  reg `var' i.opcat $polscov, cluster(fid)
	estadd local mfe "$\checkmark$"

	// FE estimates
eststo baselinefe`var':  xtreg `var' ib2.opcat $fecov, vce(cluster fid) fe
	estadd local fe "$\checkmark$"
}
  
*** Table I1: Pooled OLS estimator
esttab baselinepolsstaffratio baselinepolsnurseratio baselinepolsstaffedu baselinepolsplan baselinepolssatisfaction using "$resultspath/tablei1_baselinepols.tex", replace p starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) drop(*mc) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") compress nogap  noobs stats(mfe r2 N N_clust, fmt(0 2 0 0) labels("Municipal FE" "r$^{2}$" "Obs." "Facilities")) booktabs fragment t(3) b(3) label nobase



*** Table I2: within estimator (facility FE)
esttab baselinefestaffratio baselinefenurseratio baselinefestaffedu baselinefeplan baselinefesatisfaction using "$resultspath/tablei1_baselinefe.tex", replace p starlevels(\dag 0.1 * 0.05 ** 0.01 *** 0.001) mtitle("Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction" "Staff density" "Nurse density" "Staff education" "Plan" "Satisfaction") compress nogap  noobs stats(fe r2 N N_clust, fmt(0 2 0 0) labels("Facility FE" "r$^{2}$" "Obs." "Facilities")) booktabs fragment t(3) b(3) label nobase
********************************************************************************





*** J. Model specification tests ***********************************************

*** Table J1: Breuch-Pagan LM test of unobserved facility-specific effects
use $dataset, clear
drop if opcat == 0
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
di "`var': Breusch and Pagan Lagrangian multiplier test for random effects"
qui xtreg `var' i.opcat $polscov, re cluster(fid)
xttest0
est store bplm`var'
local bplmp`var' =  r(p)
local bplmchi2`var' = r(lm)
}
	// table
matrix bplm = (`bplmchi2staffratio',`bplmpstaffratio' \ `bplmchi2nurseratio',`bplmpnurseratio'  \ `bplmchi2staffedu',`bplmpstaffedu' \ `bplmchi2plan',`bplmpplan' \ `bplmchi2satisfaction',`bplmpsatisfaction')
matrix colnames bplm = "BP-LM-statistic" "p"  
matrix rownames bplm = "Staff density" "Nurse density"  "Staff education" "Plan" "Satisfaction"
matlist bplm
esttab matrix(bplm, fmt(2 3)) using "$resultspath/tablej1_bplm.tex", replace booktabs fragment nomtitles



*** Table J2: Hausman test comparing RE and FE estimates
use $dataset, clear
drop if opcat == 0
foreach var of varlist staffratio nurseratio staffedu plan satisfaction {
di "`var': Hausman test of Fixed effects"
qui xtreg `var' i.opcat $polscov, re
est store re`var'
qui xtreg `var' i.opcat $fecov, fe
est store fe`var'
hausman fe`var' re`var', sigmamore
local hausmanp`var' =  r(p)
local hausmanchi2`var' = r(chi2)
di "hausmanp`var'"
}
	// table
matrix hausman = (`hausmanchi2staffratio',`hausmanpstaffratio' \ `hausmanchi2nurseratio',`hausmanpnurseratio'  \ `hausmanchi2staffedu',`hausmanpstaffedu' \ `hausmanchi2plan',`hausmanpplan' \ `hausmanchi2satisfaction',`hausmanpsatisfaction')
matrix colnames hausman = "Chi$^{2}$" "p"  
matrix rownames hausman = "Staff density" "Nurse density"  "Staff education" "Plan" "Satisfaction"
matlist hausman
esttab matrix(hausman, fmt(2 3)) using "$resultspath/tablej2_hausman.tex", replace booktabs fragment nomtitles
********************************************************************************





